Non-equilibrium dynamics of a driven Bose-Einstein condensate at finite temperature 
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While the Gross-Pitaevskii equation is well-established as the canonical dynamical description of atomic 
Bose-Einstein condensates at zero-temperature, describing the dynamics of Bose-Einstein condensates at fi- 
nite temperatures remains a difficult theoretical problem, particularly when considering low-temperature, non- 
equilibrium systems in which depletion of the condensate occurs dynamically as a result of external driving. 
The BEC analog of the quantum ^-kicked rotor is the prototypical example of such a system. In this paper, 
we describe a fully time-dependent numerical method to propagate the equations of motion of a second-order, 
number-conserving description; these equations describe the coupled dynamics of the condensate and non- 
condensate fractions in a self-consistent manner, and correctly capture the phonon-like nature of excitations at 
low temperature, making them ideal for the study of low-temperature, non-equilibrium, driven systems. We 
use this numerical method to systematically explore the finite-temperature dynamics of the (5-kicked-rotor-BEC. 
We demonstrate that several qualitative features of this system at zero temperature are well-preserved at finite 
temperatures, and predict a finite-temperature shift of resonance frequencies which could be verified by future 
experiments. 
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I. INTRODUCTION 

Even at zero temperature, typical atomic Bose-Einstein 
condensates (BECs) contain a small non-condensate fraction 
due to inter-atomic interactions. In real experiments, which 
necessarily take place at finite temperatures and may involve 
dynamical depletion of the condensate [1] — for instance due 
to non-equilibrium dynamics induced by changes in applied 
external fields, a non-negligible non-condensate fraction often 
arises. A full understanding of the dynamics of such systems 
requires a fully-dynamical, finite-temperature theoretical de- 
scription which goes beyond the mean-field, zero-temperature 
Gross-Pitaevskii equation (GPE). Due to the complexity of 
such descriptions, the non-equilibrium dynamics of atomic 
BECs in the presence of a significant non-condensate fraction 
(whether thermal or dynamical in origin) remains a largely 
open problem [2, 3]. 

In systems where the non-condensate fraction is primar- 
ily thermal in origin, a variety of theoretical descriptions 
have been developed and successfully applied to atomic 
BEC dynamics. These include symmetry-breaking descrip- 
tions, which are based on a perturbative expansion about a 
mean field; for example the Hartree-Fock-Bogoliubov-Popov 
(HFBP) description [4-7], and the Zaremba-Nikuni-Griffin 
(ZNG) description [8, 9] which has been used to study ele- 
mentary excitations [10, 11], vortices and vortex nucleation 
[12, 13], and dark solitons in trapped condensates [14]. Other 
successful descriptions have been obtained in the context of 
c-field methods, which describe the highly occupied modes of 
the system as a classical field. These descriptions include the 
projected Gross-Pitaevskii equation (PGPE) [15-17], which 
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has proved successful in studies of the critical temperature 
Tc [18], and vortex formation [19]; the truncated Wigner 
PGPE, [20] which has been applied to diverse problems such 
as vortex lattice formation [21], colliding [22], reflecting [23] 
and collapsing [24, 25] condensates, collisions of bright soli- 
tary waves [26], and BEC interferometry [27, 28]; and the 
stochastic projected GPE (SPGPE) [29-32] and stochastic 
GPE (SGPE) [2, 33-37] which have been used to study the 
condensation phase transition [32, 38], and vortex dynamics 
and superfluid turbulence [32, 39, 40]. 

Alternative systems, in which dynamical depletion of 
the condensate is most important, generally involve a low- 
temperature BEC being driven by an applied external field. 
Such systems include atomic BEC analogues of generic quan- 
tum chaotic systems; e.g., the kicked accelerator [41, 42], 
kicked harmonic oscillator [43-45], and kicked rotor [46- 
51]. These systems offer an excellent test bed for explor- 
ing generic issues of quantum chaos [48, 52], quantum su- 
perposition [53], quantum resonances [43, 49, 50], dynami- 
cal instability and dynamical depletion [45, 46, 49, 50], and 
even entropy, thermalization and integrability [54-56]. In or- 
der to comprehensively describe such systems one must cor- 
rectly capture the interplay between driving, condensate, and 
low-lying non-condensate excitations. Consequently, using a 
theoretical description which self-consistently includes such 
effects is of paramount importance. Unfortunately, the finite- 
temperature descriptions mentioned above are not suited to 
such systems: the HFBP and ZNG descriptions incompletely 
model the phonon character of elementary excitations at low 
temperatures, the PGPE is restricted to a high-temperature 
regime, and truncated Wigner PGPE is prone to spurious ther- 
malization [20] which renders it accurate only for short times 
and close-to-equilibiium dynamics [57]. At low temperature 
the SPGPE reduces to the truncated Wigner PGPE, and suffers 
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from the same problem. 

In contrast, number-conserving descriptions [58-62], in 
which the system is partitioned into manifestly orthogo- 
nal condensate and non-condensate parts using the Penrose- 
Onsager criterion [63], are, in principle, ideally suited to this 
regime. However, until recently only the first-order number- 
conserving description of Gardiner [59] and Castin and Dum 
[58] offered a description which was numerically tractable 
in fully time-dependent form. Application of this first-order 
number-conserving description to the 5-kicked-rotor-BEC 
[46, 47, 49] and (J-kicked-harmonic-oscillator-BEC [45, 64] 
revealed a general tendency for rapid, unbounded growth. Un- 
fortunately, as the equation describing the condensate in this 
first-order treatment is simply the GPE, this rapid, unbounded 
growth can be viewed as a consequence of linearized instabili- 
ties present in the GPE, which would be removed by a higher- 
order description containing a self-consistent back-action of 
the non-condensate on the condensate [1]. Thus, the first- 
order description is of limited use in describing the long-time 
dynamics of low- temperature, driven BEC systems. 

To incorporate such a self-consistent back-action, a second- 
order number-conserving description is required. Such a de- 
scription was presented by Gardiner and Morgan in Ref. [65], 
and was previously used highly successfully by Morgan [66- 
68] to calculate, using a linear response treatment, the excita- 
tion frequencies of a BEC at finite temperature as measured in 
experiments at JILA [69, 70] and MIT [71, 72]. More recently 
[51], the first fully-dynamical implementation of the second- 
order equations of motion in this description was applied to 
an initially zero-temperature 5-kicked-rotor-BEC. This appli- 
cation revealed that the second-order description does indeed 
contain a self-consistent back-action which damps out the 
rapid initial growth seen in the first-order description. Con- 
sequently, this second-order description provides an excellent 
tool for studying low-temperature, driven BEC systems. 

In this paper we present, in detail, a numerical method 
for evolving the second-order number-conserving equations 
of motion at finite temperatures, which enables systematic ex- 
ploration of non-equilibrium BEC dynamics at finite temper- 
atures. A restricted implementation of this method, limited 
to zero-temperature initial conditions, was used to obtain the 
results in Ref. [51], but not described in detail in that work. 
Here, we use the finite-temperature numerical method to sys- 
tematically explore the non-equilibrium dynamics of the 6- 
kicked-rotor-BEC at finite temperature, extending the results 
of Ref. [51]. Remarkably, we observe that many of the fea- 
tures observed at zero temperature in Refs. [50, 51] are qual- 
itatively preserved at temperatures sufficient to cause signif- 
icant condensate depletion. In particular, our exploration of 
the effects of initial temperature predicts a finite-temperature 
shift in the system's resonance frequencies which could be ex- 
perimentally measured and verified. 

We begin by introducing, in Section II, the second-order 
number-conserving description of the dynamics and present 
the second-order equations of motion in their most general 
form. We give a detailed discussion of the self-consistent 
properties and regime of validity of the resulting description 
at both zero and finite temperatures. In particular, we dis- 



cuss potential problems finding equilibrium initial conditions 
at higher temperatures. In Section III we develop a numerical 
method to evolve initial states according to the second-order 
equations of motion. This method explicitly includes the non- 
linear, non-local terms coupling the condensate and the quasi- 
particle modes. It is these terms which maintain orthogonal- 
ity between the condensate and non-condensate, and which 
are particularly difficult to deal with in the context of general 
numerical integration schemes. In Section IV we apply the 
method to systematically explore the effects of finite initial 
temperatures in the 5-kicked-rotor-BEC. We also extend the 
analysis of this system given in Ref. [51] by further exploring 
the contrast between first- and second-order descriptions, and 
analyzing the evolution of the single-particle von Neumann 
entropy. We predict a shift of resonance frequencies at finite 
temperature which could feasibly be detected in experiments. 
Section V comprises the conclusions. 



II. SECOND-ORDER, NUMBER-CONSERVING 
THEORETICAL DESCRIPTION 



A. Motivation 

We consider a system of bosonic atoms of mass m, con- 
fined by an external potential V{r, t), and interacting via pair- 
wise i-wave contact interactions. The Hamiltonian for such a 
system is given by 



^(r), (1) 



where "^{r) and ^P'Xr) are second-quantized field operators 
obeying standard equal-time bosonic commutation relations. 
Here, the single-particle Hamiltonian is 



H,^{r,t)^-—^' + V{r,t), 



where 



(2) 



(3) 



for i-wave scattering length a,. 

Computing the full many-body dynamics of such a system 
directly from the Hamiltonian Eq. (1) is, in general, an in- 
tractable problem. However, in the case where is large and 
the majority of the system is Bose-Einstein condensed, one 
can obtain an approximation to the full many-body dynamics 
via a perturbative description, in which the non-condensate 
fraction constitutes the small parameter. To develop such a 
description, we adopt the definition of Bose-Einstein conden- 
sation — applicable to a finite-size and interacting system as 
described by Eq. (1) — given by Penrose and Onsager [63]. 
In this definition the condensate mode, 0c(r, f) is identified 
as a single macroscopically-occupied eigenstate of the single- 
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(4) 



particle density matrix' 

p(r,r',f) = (^plXrO^PCr)) . 
Consequently, the condensate mode </ic(r, f) satisfies 

drpir, r', O^r', t) = A^,(0<^c(r, t) , (5) 



where the condensate mode occupation, Ndt), is taken to be 
much greater than all other single-particle mode occupations 
[given by the other eigenvalues of p(r, r', ?)] at any time. 

Using this definition, one can explicitly partition the field 
operator into a condensate part, <pc{r, t), and a non-condensate 
part, S^{y, t): 



T(r) = fl,(f)^*c(r,f) + 5'^'(r,f)• 



(6) 



Here, the condensate creation operator adt) creates an atom 
in the condensate mode 0( (r, t). The hermiticity of the single- 
particle density matrix ensures that the condensate mode 
0c(r, f) is explicitly orthogonal to the non-condensate mode 
at all times. This partition is number-conserving, in the sense 
that the system remains in a state of fixed total atom num- 
ber This is in contrast to a symmetry-breaking perturbative 
description, which is based on a partition of the field opera- 
tor into a finite expectation value, ^^(r, f) = ^4'(r)^, and an 
operator-valued fluctuation §(r, t). 



*P(r) = T(r, f) + 6{r, t) , 



(7) 



which breaks the overall U(\) phase-symmetry of the Hamil- 
tonian. Note also that in such a symmetry-breaking descrip- 
tion, ^^(r, t) and 6{r, t) are only explicitly orthogonal in the 
homogeneous case. This leads to considerable difficulty in 
interpreting beyond-mean-field non-equilibrium dynamics in- 
volving both the condensate and non-condensate within a 
symmetry-breaking description; the avoidance of this issue 
constitutes a key motivation for the number-conserving ap- 
proach. 

A number-conserving description of the system's dynamics 
is obtained by treating the non-condensate part S^{r, t) pertur- 
batively, by means of an expansion in terms of a suitable fluc- 
tuation operator [58, 59, 65]. Neglecting the non-condensate 
completely yields a zeroth-order description in the form of 
the Gross-Pitaevskii equation (GPE); this takes the form [we 
explicitly denote the condensate mode in this zeroth-order de- 
scription by 0f ''(r, t) for clarity in later sections] 



.(0) 



(r,f) 



dt 



= [//,p(r, t) + U\cl>f\r, tf - <pf(x, t) , (8) 



where U - UqNc, and is given by 



<ir(^*"' (r,f) 



H,pir,t) + U\<pf>ir,tr-ih 



(9) 



' Note that in Eq. (4) the brackets (■ ■ ■ > denote an expectation value with 
respect to the full many-particle density matrix; at finite-temperature this 
involves a thermal, as well as a quantum, average. 



As in Refs. [51, 59, 65], /Iq^ is defined in such a way as to be 
generally time-dependent. However, as A^^ is explicitly real, 
the dynamics resulting from its time evolution consist only 
of an irrelevant global phase. When considering, as we do 
in this paper, the evolution of a stationary, equilibrium initial 
condition subject to a time-dependent perturbation, it is most 
convenient to work with a time-independent GPE eigenvalue 
/l!!'' given by 



Af = J dr^f^'ir, 0) f//,p(r, 0) + f7|0f >(r, O)!^! ^^(r, 0) , 

(10) 

where ^['''(r, 0) represents the t - Q stationary, equilibrium 
state of the GPE. To work with such an eigenvalue in the GPE 
description, one simply makes the replacement /ij,'* — > /l™^ in 
Eq. (8). 

In contrast to the GPE, the second-order number- 
conserving description we use in this paper provides a descrip- 
tion of both the condensate and the non-condensate, and con- 
sists of mutually-coupled equations for both, which we outline 
in the following section. 



B. Equations of motion 

Conducting a self-consistent, second-order, number- 
conserving expansion — as detailed in Ref. [65] — leads to a 
number-conserving generalized GPE (GGPE) for the dynam- 
ics of the condensate mode <^c(r, t): 



iti 



dt 



[//sp(r)-4"]0.(r) 



H- U 



h(r, r) 



N^^'^'^ ^' N, 



PAr) + U<p:(r) 



m(r, r) 

Nr 



-U j dr'\Mr'f(^ 



|0.(r')|^(%^0.(r').C(r')^]. dD 



Here we have introduced the convention, adopted gener- 
ally hereafter, of omitting explicit time-dependences to aid 
clarity. The dynamics of the non-condensate enter the 
GGPE [Eq. (11)] through the normal [n(r, r')] and anomalous 
[m(r, r')] averages. These are respectively defined by 



and 



where 



n(r,r') = (A"'"(r')A(r)) , 
m(r,r') = (A(r')A(r)) , 



A(r) 



JNr 



(12) 



(13) 



(14) 



is the number-conserving fluctuation operator in which a per- 
turbative expansion has been conducted. Note that the cor- 
responding small parameter is -\/N,(t)/Nc{t), where A^,(f) = 
- Nc{t) is the occupation of the non-condensate. We have 
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also introduced the complex, time-dependent GGPE "eigen- 
value" /Ij', which is given by 



9 m(r, r) 



j drcf>:(r)iM, 



.sp(r) + U 

h(r, r) 



+2- 



-/;!-^0,(r). (15) 



It is essential to emphasize that this "eigenvalue" is in general 
time-dependent, and that Eq. (15) should be taken to be eval- 
uated using the values of 0c(r), Nc, etc. at a particular time. 
As we demonstrate in the subsequent section, this time depen- 
dence of ^2 is crucial in order to obtain self-consistent num- 
ber dynamics. However, as in the case of the GPE [Eqs. (8) 
and (10)] it is considerably more convenient to work with a 
real, time-independent eigenvalue, and in Section III we will 
reformulate the second-order equations of motion in order to 
replace A^^ with such a quantity, A^^\ 

The GGPE must be coupled to a consistent set of equa- 
tions of motion for the non-condensate; at second-order these 
consist of the number-conserving modified Bogoliubov-de 
Gennes equations (MBdGE); 



5f I A'(r) 



, Z{r,r') M(r,r') \ A(r') 
'^'^ '-M*(r,r') -r(r,r') UV)'' ^^^^ 



where 

Ur, r') = 6{r - r') [H,^{r') + UWr')f - 

-H j dr"Q(r,r")U\cf>,(r"fQ(r",r'), (17) 

and 

yVl(r,r') = J dr"Qir,r")U(f>ir"fQ*ir",r'). (18) 

The modification in these equations — with respect to the or- 
dinary BdG equations obtained in a symmetry-breaking de- 
scription — is the appearance of the projector (3(r, r'), defined 
by 



Q(r,r') = 5(r'-r)-0,(r)C(r'), 



(19) 



which explicitly enforces the orthogonality of the condensate 
and non-condensate, and of the "GPE eigenvalue" Aq\ This 
"GPE eigenvalue" is obtained by substituting the GGPE con- 
densate mode 0£.(r, f) into Eq. (9) in place of the GPE conden- 



sate mode 0f '*(r, f). That is. 



drfjr, t) 



H,^{r,t) + U\<t>c{r,tt -in 



dt 



(20) 



where we have ressurected time-arguments for clarity. As in 
the zeroth-order, GPE, case — where it was possible, and con- 
venient for the case of a stationary, equilibrium initial state, to 
replace the time-dependent GPE eigenvalue /Iq' with a time- 
independent GPE eigenvalue /l^'* — is an explicily real 
quantity. It will thus later be convenient, for the case of a sta- 
tionary, equilibrium initial state, to replace the time-dependent 
"GPE eigenvalue" with the time-independent "GPE eigen- 
value" Ig"' given by 



1(0) 



J drf^ir, 0) [i/,p(r, 0) + U\Ur, 0)|'] <^.(r, 0) , (21) 



where we have again resurrected the (zero) time-arguments 



for clarity. Note, however, that the replacement of /ig' with 



A^q'' must be accomplished consistently with the replacement 
of /Ij' with /Ij*^' as outlined in Section III. 

The generalized Gross-Pitaevskii equation [Eq. (11)] and 
the modified Bogoliubov-de Gennes equations [Eq. (16)] 
complete the fully dynamical, second-order, number- 
conserving description obtained by Gardiner and Morgan in 
Ref [65] and previously used within a linear response treat- 
ment by Morgan [66-68]. In this paper, we develop these 
equations into a form where their fully dynamical time evolu- 
tion can be realized numerically, as in Ref. [51]. 



C. Discussion 

Before we outline our method for the simultaneous numer- 
ical solution of Eqs. (11) and (16), a few comments regarding 
the properties of the second-order equations of motion and 
their regime of validity are in order. Firstly we note that, in 
general, the diagonal part of the anomalous average, fii{r, r), 
is ultraviolet-divergent and must be appropriately renormal- 
ized; this procedure is described in detail in Refs. [65] and 
[62]. One exception is in the case of quasi-one-dimensional 
systems, such as the one we consider in Section IV, where 
renormalization is not necessary. 

Secondly, the number-conserving equations of motion used 
here are derived by expanding the total Hamiltonian in powers 
of the number-conserving fluctuation operator A(r) [Eq. (14)]. 
This operator is advantageous for three primary reasons: (1) it 
scales proportionally to the number of non-condensate atoms, 
which we wish to treat as a small parameter; (2) it avoids the 
need to expand inverse-square-root number-operators when 
expanding the Hamiltonian, which is particularly convenient; 
and (3) it is a well-defined fluctuation operator in the sense 
that ^A(r, f)^ = 0. It should be noted that the commutation 
relations of A(r) 



[A(r),A^(r')l=^ar,r') 



1 



(22) 



6A>\r')S^{r) , 



give rise to quasiparticles which are, in principle, only approx- 
imately bosonic. However, when restricted to the (quadratic) 



order of approximation we consider, the corresponding quasi- 
particles are indeed exactly bosonic. Hence, the favorable 
properties of A(r) mentioned above make it a preferable 
choice for developing the second-order equations of motion 
[62, 65-68, 73]. 

Another potential issue is that the appearance of higher- 
than-second-order terms in the equations of motion has been 
prevented by working within a consistent Gaussian approxi- 
mation [65]; that is, all quadratic products of operators are 
required to take the form of pair averages. This constitutes 
a Gaussian approximation in that all higher-order moments 
of the fluctuation distribution are assumed to be describable 
in terms of (A"'"(r)A(r)) and (A(r)A(r)) [74]. To enforce 
this requirement, cubic products of fluctuation operators have 
been factorized into products of single fluctuation operators 
multiplied by pair averages. This approximation is identi- 
cal to the Hartree-Fock factorization of similar cubic terms 
in symmetry-breaking treatments, and corresponds to making 
the replacements 

A"'"(r)A(r')A(r') ^ 2 (A"'"(r)A(r')) A(r') 

+ (A(r')A(r'))Af(r), (23) 

etc. As discussed in detail in Refs. [2, 62, 65], this proce- 
dure potentially leads to an inconsistent treatment of interac- 
tions. However, as was explicitly demonstrated by Morgan 
[62], the primary problem with Hartree-Fock factorization of 
cubic operator products is that it can lead to the omission 
of cubic terms which are actually larger than quartic terms 
which are subsequently retained. In the second-order number- 
conserving description presented here this inconsistency does 
not arise, as all quartic terms are neglected. 

A key property of the second-order number-conserving 
equations of motion is their number-self-consistency: in con- 



(a) 



trast to aI?^ and X 



3(0) 



which can both be considered as low- 



Af is a 



'0 

order approximations to the chemical potential 
complex eigenvalue. The meaning of the imaginary part of 
/Ij^ can be understood by considering the (implicit) time de- 
pendence of Nc, which is given (to quadratic order) by 



^ dNc ^ d 

in = in — 

dt dt 



N- Jdr (A^(r)A(r)) 




(24) 



(rfm(r,r)-ih*(r,rWrf] 



= (A''' - A 



Thus the time-dependent, imaginary part of A^^ acts to keep 
the condensate mode <pcir) normalized to unity despite the 
growth or decay of the condensate population. This illus- 
trates the presence of number-self-consistency in the dynami- 
cal coupling between the GGPE and MBdGE. In contrast, the 
first-order description [58, 59] — which consists of the same 
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FIG. 1. Schematic representation of the zeroth-, first-, and second- 
order number-conserving equations of motion. At zeroth order, (a), 
the non-condensate is ignored and the condensate mode (pf \r) is de- 
scribed by the Gross-Pitaevskii equation (GPE) [Eq. (8)]. At first 
order, (b), the GPE [Eq. (8)] is coupled to modified Bogoliubov-de 
Gennes equations (MBdGE) [Eq. (16)] for the non-condensate fluc- 
tuation operator A(r). At this order, the evolution of A(r) depends 
on the evolution of (r), but the converse is not true; this order of 
approximation can be interpreted as treating the condensate as an in- 
finite atomic reservoir. At second order, (c), the non-condensate is 
again described by the MBdGE [Eq. (16)]; however, the evolution 
of the condensate mode 0r(r) is now determined by the generalized 
Gross-Pitaevskii equation (GGPE) [Eq. (11)]. This pairing of equa- 
tions produces fully self-consistent number dynamics. 



MBdGE [Eq. (16)] coupled to the ordinary GPE [Eq. (8)] — 
is not in general number-self-consistent, since the conden- 
sate population is fixed. The first-order description can only 
be considered to be number-self-consistent when viewed as 
the limit of the second-order description as ^ oo [65]. 
The zeroth-order description, consisting of the GPE [Eq. (8)] 
alone, is trivially number-self-consistent, as it ignores the 
growth and decay of the condensate altogether. The number- 
self-consistency of the zero-, first-, and second-order number- 
conserving descriptions are illustrated schematically in Fig. 1. 

A final feature of the second-order, number-conserving 
equations of motion is that the non-condensate is described 
by a MBdGE [Eq. (16)] which does not contain pair aver- 
ages of the non-condensate operators; the terms Ai{r, r') are 
in no way altered from the first-order description and the 
terms X(r, r') appearing in the MBdGE consist only of a GPE 
Hamiltonian and the "GPE eigenvalue" of Eq. (28). Impor- 
tantly, however, the GPE Hamiltonian and the "GPE eigen- 
value" are, in the second-order description, evaluated in terms 
of the second-order GGPE wavefunction <;*c(r). This leads to 
an apparent inconsistency in that the spinors [0c(r), 0]^ and 
[Q,<pl(r)]^ are no longer exact, zero-energy solutions of the 
MBdGE, as they would be for the GPE wavefunction (p^^K 
This has the unfortunate side-effect of making the identifica- 
tion of self-consistent initial conditions difficult at high tem- 
peratures, particularly in inhomogeneous systems. However, 
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this problem can be viewed purely as a consequence of ap- 
plying the theory outside of its regime of validity, as argued 
by Morgan in his study of condensate excitations at finite- 
temperature [66-68]; in this work he demonstrated that the 
second-order description remains self-consistent at high tem- 
peratures (approaching A^, ~ Nc) provided one restricts one- 
self to a linear response treatment, in which a self-consistent 
equilibrium solution is not necessary. A fully dynamical 
treatment, which requires a self-consistent equilibrium initial 
condition, is thus restricted to lower temperatures (such that 
Nf < Nc) than a linear response treatment. 



III. FULLY TIME-DEPENDENT NUMERICAL 
IMPLEMENTATION 

A. Reformulation of equations of motion 

1. Elimination of complex, time-dependent "eigenvalue" 

In this Section we develop a numerical method for evolving 
the combined GGPE and MBdGE system of Eqs. (11) and 
(16). In order to do so, it is of great convenience to ehminate 
the imaginary part of the GGPE "eigenvalue" A^^ [Eq- (15)]; 
this can be done by describing the condensate using a mode 
function normalized to the condensate population, Nc- Hence, 
we define 



and the adjusted GGPE eigenvalue. 



i/r(r)= y/Nc<pc(r), 



(25) 



in terms of which the GGPE can be re-expressed [using 
Eq. (24) for the number evolution] as 



= K(r) + t/ol«A(r)|' - ] ^(r) 



h(r, r) - 



\>Hrf 



+ f/o«(r,r)i/r(r)- 



N 



+ Uofhir, r)(A*(r) ■ 



and the MBdGE as 



N 



"J dr' i/r(r' 
^ J dr'rir' 



tA(r) 

)|(A(r')l'«(r,r') 



Mir'fmry), (26) 



iti-Mr) = [//,p(r) + f/ol<A(r)l' - io*°'] A(r) 

+ t/o|iA(r)|'A(r) - ^ J dr' r{r'mr'fk{r')m 

^ J fl'r'iA(r')|«A(r')l'A"'"(r')'A(r). (27) 



+ f/o-A(r)'A"'"(r) 



Here, we have chosen to write the equations in a form which 
emphasizes the significant structural analogies between the 
GGPE and the MBdGE. The eigenvalues appearing in these 
equations are the "GPE eigenvalue", 

= ^ J fl^r ^*(r) [i/,p(r) + Um^f] <A(r) , (28) 



,(0) 



^-Af.'^Jdrrir) 



2«(r,r)-— |iA(r)P 



(A(r) 



2N 



- [iA(r)*^OT(r, r) + i^(rfm(r, r)] . (29) 



Both of these quantities (Eq. (28) and Eq. (29)) are now ex- 
plicitly real, and should be evaluated explicitly in terms of the 
stationary, equilibrium initial condition of the GGPE at f = 0. 



2. Quasiparticle decomposition 

In numerical studies of non-equilibrium dynamics at fi- 
nite temperature one generally wishes to begin from a finite- 
temperature equilibrium condition and explore the dynamics 
resulting from, e.g., driving or an abrupt quench event in 
a fully-time dependent way-. In the second-order number- 
conserving description such a thermal and dynamical equilib- 
rium state corresponds to a self-consistent, stationary solution 
of the equations of motion in which the elementary quasipar- 
ticle excitations of the system are populated according to the 
appropriate thermal Bose distribution. 

The appropriate quasiparticle basis in the second-order 
number-conserving description is the basis of Bogoliubov 
quasiparticle modes which diagonalize the stationary MBdGE 
[65]: 

/ X(r,r') M(r,r') \ \^ ^ I Uk(r)\ ,,,.\ 

-Z^4"f«)^"''*^'''^'"'^'''^^ ■ ^^^^ 



In terms of the quasiparticle mode functions Mt(r) and vt(r), 
the fluctuation operator A(r) can be expanded as 

where hi and bk are quasiparticle creation and annihilation 
operators which, at this order, are assumed to have bosonic 
commutation relations (see Section II C). Assuming all time- 
dependence to reside in the quasiparticle mode functions Mjt(r) 
and vt(r), with the quasiparticle creation and annihilation op- 
erators and b^ time-independent, the MBdGE [Eq. (27)] 
take the form 

in-ukir) = [//,p(r) + f/ol<A(r)l' - A^] Uk{r) 

+ Umrtut{r)-^ j dr' r(r'M(r')fuk(r')m 
+ Uom\(r)-Y Jdr'il,(r'mr')\\(r'),p(r), (32) 



^ One could in principle start from any non-equilibrium initial condition, but 
we consider only the equilibrium case here. 
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and 



injvl(r) = [//,p(r) + f/„|-A(r)|' - ~Af>] vl{r) 

+ Uomr)fvl(r) - ^ J dr' r(r'M(r')fvl(r')i^(r) 

+ UoHrful(r)-^ Jdr'i/,(r'mr'ful(r'Mr). (33) 

At initial thermal equilibrium, the quasiparticle creation 
and annihilation operators have the following pair averages: 



{bib,) = 6uNk , 
{hb)^{blb\)^Q, 



where A^^ is the Bose-Einstein factor [62, 65, 67] 

-1 



exp 



1 



(34) 
(35) 



(36) 



The term ^ - A^^^ represents a finite-size correction, given by 
[62, 67, 75] 



(37) 



This leads to quasiparticle expressions for the non-condensate 
normal and anomalous averages [Eqs. (12) and (13)] 

oo oo 

«(r, r') = 2 NkUk(r)ul(r') + ^(A^t + l)v*(r)vi(r') , (38) 

k=\ k=l 

CO OO 

m(r, r') = ^ A^,«i(r)y*(r') + J](Nk + l)v:(r)«,(r') . (39) 



k=l 



k=\ 



We reiterate that, by choosing a scheme where all time depen- 
dence resides in the quasiparticle mode functions Uk{r) and 
Vk{r), the quasiparticle populations A^^ appearing in Eqs. (38) 
and (39) remainfixed in time. 

Using the above expressions we proceed, in the next Sec- 
tion, to re-cast the GGPE in terms of the quasiparticle mode 
functions, and re-cast the combined GGPE and MBdGE — 
now in the form of Eqs. (26), (32), and (33) — in the final 
form which we will use to conduct a simultaneous numerical 
solution. 



3. Introduction of the spinor i, 



of the MBdGE can be computed by ignoring the projector 
terms in the MBdGE throughout the evolution, and simply 
projecting the final state orthogonally to the condensate [45]. 
However, if one were to similarly ignore the projectors dur- 
ing the evolution of the second-order GGPE-MBdGE system, 
one would then have to re-orthogonalize both the condensate 
and quasiparticle modes with respect to an unknown basis at 
the end of the evolution. Consequently, the projection terms, 
and the non-local integrals they involve, must be explicitly in- 
cluded in any numerical method. In this Section we develop 
such a method, by re-casting the GGPE and MBdGE in a form 
which exploits their apparent symmetries, and allows us to in- 
clude the projection terms using a split-step technique. 

After substituting the quasiparticle expressions for the nor- 
mal and anomalous average [Eqs. (38) and (39)] into the 
GGPE [Eq. (26)], the GGPE can be recast in the form 



dMr) 

= [i/Gp(r) + B(r)] ^(r) 



dt 



+ + ^)Ak(rH(r) + 2 NkAl(r)uk(r) , (40) 

k=\ k=i 

and the MBdGE [Eqs. (32) and (33)] can be recast in the form 



in"^^ = HGp(r)uk(r) + Akir^r) , 



ifi 



dt 
dvlir) 

dt 



where 



(0) 



Hopir) = Hsp(r) + Uoimr - ~A, 
Ak(r) = Uo [vk(r)ifr(r) + Mt(r)-A*(r) - h] 

^ J flfr [v,(r)iA(r) + Uk(r)rir)] |<A(r)|' 



N, 



(41) 
(42) 



(43) 
(44) 

(45) 



and 



B(r) = Uo 



Y^Nk\uk(r)f + (Nk + l)\vk(r)f 



k=l 



-Uq- 



- A' , (46) 



where A' = A^^ - A^^K This reformulation of the problem 
allows one to write the coupled evolution of the condensate 
wavefunction and the first M quasiparticle modes as a nonlin- 
ear matrix equation in a (2M + l)-dimensional spinor space: 



The primary difficultly in simulating the coupled GGPE- 
MBdGE system numerically is the problem of orthogonaliza- 
tion: both equations contain terms which function to maintain 
orthogonality between the condensate and non-condensate. 
This is in contrast to the case of the first-order GPE-MBdGE 
system, where the GPE evolves in isolation from the MBdGE; 
this de-coupling in the first-order system means the evolution 



/;z-^(r) = r(r)^(r). 
ot 



(47) 



Here the vector ^(r) is defined by 

ar) = [^(r), v\(r), y^(r), ui(r), UM(r)f , (48) 
and the operator r(r) is defined by 



8 



(HcAr) + B(r) (A?, + l)Ai(r) (A?2 + Wr) 



r(r) 



A*(r) 
A*(r) 



Ai(r) 
A2(r) 



AM(r) 



Hcp(r) 






^/Gp(r) 







{Nm + l)AM(r) NiA;{r) A?2A*(r) 



i^Gp(r) 











Hcpir) 
//Gp(r) 



A^MA;,(r) 









HGp(r) ) 



In any actual calculation, this spinor space is rendered finite- 
dimensional by the need for a finite quasiparticle momentum 
cut-ofi" M. However, M may, in principle be arbitrarily large. 



B. Operator-splitting sclieme for time-evolution 

1. Separation of position and momentum terms 

As we have already accounted for all creation and annihila- 
tion operators through the quasiparticle decomposition, each 
entry in the matrix defining r(r) can be thought of as an op- 
erator in the first-quantized sense^^. From an analytic perspec- 
tive, this notation seems to achieve little more than 'tidying' 
— abstracting away much of the detail. Importantly, how- 
ever, all the operators which are off-diagonal in the spinor 
space [the operators Aj.(r)] are diagonal in the position rep- 
resentation, since they ultimately consist of a multiplication 
by a spatially-varying function. In contrast, all the operators 
which have off-diagonal components in the position represen- 
tation [that is, the kinetic energy operator implicitly contained 
in //sp(r) and hence in //Gp(r)] appear only on the diagonal in 
the spinor space. 

From a numerical perspective this arrangement is extremely 
useful, as it makes the evolution amenable to a split-step ap- 
proximation. This is achieved by splitting r(r) into the sum 
of a term representing the linear single-particle evolution, 
rL(r), and a term representing nonlinear parts of the evolu- 
tion, FnCf). The linear, single-particle term is defined by 

rL(r) = /®i/L(r), (50) 

where / represents the (2M -H 1) x (2M -i- 1) identity matrix, 
and we have also defined 

HlW = //sp(r) - (51) 

The nonlinear term is then defined simply by 

rN(r) = F(r)-FL(r). (52) 



^ That is, an operator which could appear on the right side of a single-particle 
Schrodinger equation. 



Note that this matrix contains diagonal entries of the form 
H^{r) and H^^ir) + B(r) where 

H^(r) = Hcpir) - H^r) ■ (53) 

Written in this form, FL(r) contains all kinetic energy terms; 
while these terms are not diagonal in the position representa- 
tion, they do all lie on the diagonal in the spinor space. Con- 
sequently, the evolution due to the kinetic energy terms can 
be computed for each mode function separately. In contrast, 
while rN(r) does contain off-diagonal elements in the spinor 
space, each element of rN(r) is diagonal in the position repre- 
sentation. Consequently, the evolution due to these terms can 
be computed straightforwardly over a discrete spatial grid. 

In the spinor space, the split-step approximation for the evo- 
lution of the system over short times 6t is given by 

^(r, t + 6t)= e-'r'(''>*'/*^(r, f) 

^ g-/rN(r)«/2Sg-iTL(r)5;//ig-irN(r)5;//i ^^^^ ^-^ ^ ^^^-^ 

where each of the three evolution operators on the right is as- 
sumed to act instantaneously, and the symmetrization reduces 
the error to the order of 6t^ [76]. Since rL(r) is diagonal in the 
spinor space, we have 



2. Analytic form for position-space evolution 

The spatial operator rN(r), is more problematic because it 
is not diagonal in the spinor space, and because it contains the 
nonlocal, nonlinear terms necessary to preserve orthogonality 
between the condensate and non-condensate. However, each 
of its entries is diagonal in the position representation, mean- 
ing we only need exponentiate FN(r) in the spinor space in 
order to obtain an operator we can evaluate and use for short- 
time propagation at a given position: such an operator — in 
distinct contrast to g^'^LCr)*'/'' which consists of an application 
of e-'fL(r)i5(/fi Q^ch mode individually — couples the mode 
functions by acting on all modes at once. The numerical ad- 
vantage of casting the evolution in terms of such an operator 
is that it almost entirely overcomes the problem of retaining 
orthogonality between the condensate and the quasiparticle 
modes, which poses severe difficulties for other schemes. In 
practice we do find that a correction for numerical round-off 



9 



error is still required in order to preserve complete orthog- 
onality over long times; this can be achieved by explicitly 
orthogonalizing the quasiparticle modes with respect to the 
condensate at the end of each time-step. However, our expe- 
rience is that applying such a correction after every time-step 
does not cause decay of the total atom number, indicating that 
the need for such a correction does indeed stem from the nu- 
merical round-off error inherent in finite-precision arithmetic. 
Indeed, the problem of loss of orthogonality between vectors 
due to finite-precision arithmetic is well-known in the context 
of, e.g., the Gram-Schmidt process [77]. 



In principle, the evolution due to the term e-'rN(r)5r/2S ^.^^ 
be obtained, for a given r, by exponentiating the matrix Tt^ir) 
numerically at each spatial grid point. However, with the nec- 
essary number of floating point operations necessary to diag- 
onalize the matrix rN(r) scaling at least as (2M + 1)^ [77], 
this leads to an algorithm which is too slow to be practicable. 
To overcome this limitation, we have computed a general an- 
alytic expression for e"'rN(r)i5'/2S f^j. arbitrary M: this reduces 
the scaling of the computational effort to (2M+ 1)^ (i.e., equiv- 
alent to a matrix-vector multiplication). This expression is ex- 
plicitly derived in Appendix A, where we also show that, at 
any specific point r (and time t) on the computational grid, 
the time evolution due to Fn can be approximated locally, af- 
ter numerical evaluation of the non-local integrals appearing 
in the functions Ait(r) over the entire position-space grid, by 
the following 2M + 1 coupled equations: 



rUr) = /sin l^j +E(r) 



X exp 



2n 



6t_ 
2h 

B(r) 



+ S(r) 



and 



rcos(r) + ^rsin(r)fl(r)-re,p(r) 
S(r) 



2(r) = Y,(2Nk + im(r) 

k=i 



, (62) 



(63) 



(64) 



Note that the quasiparticle index k in Eqs. (57) and (58) runs 
from 1 to M, and we have chosen to show time-dependences 
explicitly in Eqs. (56-58), to make the time-dependent nature 
of the expression clear. The analytic evolution scheme given 
by Eqs. (56-58), incorporated into the symmetrized split-step 
method of Eq. (54), provides a way to numerically implement 
the second-order number-conserving description of Section II. 



IV. NON-EQUILIBRIUM DYNAMICS OF THE 
(5-KICKED-ROTOR-BEC AT FINITE TEMPERATURES 



(A(r, t + 6t/2) = I r,o,(r, t) - -r,i„(r, t)B{r, t)j t//{r, t) 

-nin(r,OS(r,0, (56) 



v^(r, t + dt/2) = -A*(r, f)rsin(r, t)ifr(r, t) + T,,^(r, f)v^(r, f) 

+ r„i,(r,f)H(r,f), (57) 

Uk(r, t + 6t/2) = -Ai(r, f)7',si„(r, f)'A(r, f) + T,^p(r, t)uk(r, t) 

+ r„i,(r,f)H(r,f), (58) 



where 



M 



E(r) = ^(A^, + l)A,(r)y;(r) + A^/A;(r)«,(r) , (59) 



i=\ 



and we have defined 

' -iH^{r)6t 



7'exp(r) = exp 



7'cos(r) = cos 



2n 



B(r) 



(60) 



+ E(r) 



St 

2h 



X exp 



2n 



(61) 



A. Overview 

In this Section we apply the split-step numerical method de- 
veloped in the previous Section for evolving the second-order, 
number-conserving equations of motion — consisting of cou- 
pled GGPE and MBdGE — to the (5-kicked-rotor-BEC. This 
system consists of a quasi- ID, toroidally-trapped, repulsively- 
interacting atomic BEC driven by 5-kicks from a spatial co- 
sine potential, and serves as a prototypical example of a sys- 
tem where condensate depletion occurs dynamically as a re- 
sult of external driving. The free parameters of this system 
can be expressed in terms of dimensionless kick strength k, 
kick period Tp, and interaction strength gj-, plus the total atom 
number A^, as shown in Fig. 2 (see Section IV B for a full def- 
inition of these parameters). This system is a BEC analogue 
of the quantum 5-kicked rotor [78-82], a paradigm quantum- 
chaotic system exhibiting complex behavior as a result of the 
periodic driving [78, 79, 82]. In particular, the 5-kicked ro- 
tor exhibits quantum resonant driving frequencies, at which 
the uptake of energy from the driving potential is greatly in- 
creased, with diffusive (linear in time) expansion of the state 
in momentum space being replaced by ballistic (quadratic in 
time) expansion at resonance [78-80]. Such systems have 
previously been realized in the context of atom-optics ex- 
periments [81-84]. However, the 5-kicked-rotor-BEC offers 
a route to a range of new dynamics and phenomena, since 
the nonlinear effects of inter-atomic interactions introduce, 
on the mean-field level, the potential for true wave chaos 
[45-50, 64, 85-87]. For experimentally realistic interaction 
strengths and atom numbers, the key new dynamical features 
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FIG. 2. The (5-kicked-rotor-BEC system: A repulsively-interacting 
atomic BEC consisting of A' atoms is held in a quasi-lD, toroidally- 
shaped trap, and driven by periodic kicks of short duration t,,, which 
can be modeled as (5-impulses, from a sinusoidal driving potential. 
The free parameters of the system are the dimensionless kick strength 
K, kicking period T,,, and interaction strength gj. See Section IV B 
for a full definition of these parameters. 



of the (J-kicked-rotor-BEC, within a mean-field description, 
are the appearance of nonlinear quantum resonances (with 
no analog in the linear regime) and the appearance of sharp, 
asymmetric cut-offs in the resonance profiles as a function of 
driving period. Both these features were identified in Ref. [50] 
using the zeroth-order, GPE description. 

However, the beyond-mean-field effects of condensate de- 
pletion and condensate - non-condensate interactions can be 
expected to play an increasing role for longer times in driven 
systems such as the (5-kicked-rotor-BEC. As previously stated, 
the explicit dynamical separation of condensate and non- 
condensate inherent in number-conserving descriptions makes 
them ideal for the study of such beyond-mean-field effects 
in driven BEC systems. Unfortunately, the lack of a self- 
consistent back-action of the non-condensate on the conden- 
sate in the first-order number-conserving description leads to 
unbounded growth of the non-condensate, and unphysical re- 
sults at long times when close to resonance [46, 47, 49]. The 
same effect has also been observed in the similar (5-kicked- 
harmonic-oscillator-BEC [45, 64]. In contrast to the first- 
order description, the presence of a self-consistent back-action 
in the second-order number-conserving description makes it 
ideal for the description of such systems. 

In Ref. [51] a restricted implementation of the numerical 
method of Section III, limited to the study of zero-temperature 
initial conditions, was used to demonstrate explicitly that the 
self-consistent back-action of the non-condensate damps out 
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FIG. 3. Zero-temperature evolution of the (5-kicked-rotor-BEC 
close to a nonlinear resonance explored in [50] (dimensionless kick 
strength k = 0.5, kick period Tp = 6.12, and interaction strength 
gr = 2.5 X lO"'*, and total atom number N = 10^ — see Sec- 
tion IV B). Beginning with a zero-temperature initial condition, three 
descriptions are used to determine the dynamics: the second-order 
number-conserving description [GGPE/MBdGE, Eqs. (11) and (16)], 
first-order number-conserving description [GPE/MBdGE, Eqs. (8) 
and (16)], and a first-order number-conserving description with ad- 
hoc renormalization. In each case we show the evolution of the 
condensate fraction = N^/N and the non-condensate fraction 
n, = N,/N. Note that, in the first-order description without renor- 
malization, «£. = 1 for all times. 



initially rapid growth in the non-condensate, which would 
continue unbounded in the first-order description. We il- 
lustrate this further in Fig. 3, which shows the number- 
evolution obtained using the first- and second-order number- 
conserving descriptions to evolve the dynamics of an initially 
zero-temperature (5-kicked-rotor-BEC, for parameters close to 
a nonlinear resonance. In addition to the evolution of the 
first- and second-order descriptions, which were originally ob- 
tained in Ref. [51], we show the evolution of the same initial 
conditions under "renormalized" first-order equations of mo- 
tion; these are identical to the first-order equations of motion, 
but with an added renormalization of the condensate, applied 
by explicitly calculating A^, and setting Nc - N - N, imme- 
diately after each time step. While it does prevent the occur- 
rence of unphysical growth of the total particle number, this 
ad-hoc renormalization procedure fails to damp out the rapid 
growth of the non-condensate in a similar way to the second- 
order description. This illustrates that the self-consistent na- 
ture of the back-action of the non-condensate in the second- 
order description is key to preventing unphysical growth of 
the non-condensate. In the following Sections we use the nu- 
merical method of Section III to explore the dynamics of the 
5-kicked-rotor-BEC for finite-temperature initial conditions. 
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B. Theoretical model and initial conditions 

As shown in Fig. 2, we consider a system described 
by Hamiltonian [Eq. (1)] with a toroidal external potential 
Vt(p,z) = ma?[(p - Rf + z^]/2 [Fig. 2(a)]. Potentials 
similar to this have been experimentally realized in, e.g., 
Refs. [88, 89]. We assume a sufficiently strong effective trap 
frequency oj that the harmonic oscillator length in this direc- 
tion, flr = s/H/ma) R. Under this assumption, a quasi- ID 
approximation is justified and yields the ID model Hamilto- 
nian [90] 



1 



+ v{e,t) + ^4'^(6')4'(6') 



no), 



(65) 

where we have chosen to use length units of R, time units of 
mR^/H, and we have introduced the dimensionless interaction 
strength gj = losR/af.. We restrict our analysis to the case 
where the system size, 2jtR, is much smaller than the phase 
coherence length - 2nR exp( yjnNJlg^)! ^AngjNc', in this 
case the ground state of the system can contain a true homo- 
geneous Bose-Einstein condensate [91]. We do not consider 
the alternative case of a quasicondensate (Z^ < 2nR) [92, 93]. 

As in previous studies of the (5-kicked-rotor-BEC [46, 47, 
49, 50], we model the driving potential as a train of (5-kicks 



V{0, f) = KCO?.{6) ^ S(t - riTp) . 



(66) 



n=0 



Here the dimensionless kicking period Tp is defined in terms 

(real) 



of the real kick period Tp as 



HT^r'^/mR^, and the 
dimensionless kick strength k is dependent on the real strength 
of the driving potential, and its duration tp. Such a driving 
potential may be approximated in experiment using, e.g., short 
pulses of off-resonant laser light [80-82]. Further details of 
the physical aspects of the system can be found in Refs. [50, 
51]. 

The equations of motion for the 5-kicked-rotor-BEC are 
identical to those of Section 11 B with the replacement r — > 0. 
The dynamical and thermal equilibrium state of the system 
consists of a uniform condensate mode, t//(ff) = yfNJln, ac- 
companied by a thermal population of Bogoliubov quasipar- 
ticle excitations. Both at T = and for T > Q, the stationary 
initial Bogoliubov modes associated with the uniform conden- 
sate are given by [51] 



where 



Uk(ff) 

Vkiff) 



Ck + q 

Ci - Ct 



JkB 



V2^' 



1/4 



and, for the uniform initial condensate, Ap** = gTNcl2n. 
is important to note that in this Section we have replaced 
the completely generic quasiparticle index k - 1,2, .. . used 
in Section 111 with the quasiparticle momentum index k — 



(67) 



(68) 



It 



. . . , -2, -1,1,2,..., as this considerably simplifies the nota- 
tion. Consequently, in our numerical treatment we must also 
replace the generic quasiparticle cut-off momentum M with an 
explicit maximum quasiparticle momentum index fcmax- The 
expressions appearing in Section 111 can be reformulated in 
this new notation simply by making the replacements 



M 



k— 1 k— /Tniax 

M 2k^„^ , 



(69) 
(70) 



where 2' indicates a summation omitting the term k = Q 
(which here corresponds to the condensate mode). 

The energy of a quasiparticle of momentum k is given by 



ek 



4 



(71) 



Hence, the non-condensate population can be determined 
from 

A^, = doY, [Nk\uk(9t + (Nk + l)\vk(0t] 

''■max 

k—k 

where the populations A^^. are obtained from the energies ek us- 
ing the appropriate thermal Bose distribution [Eq. (36)]. Note 
that in our system of dimensionless variables, the dimension- 
less temperature T used to obtain populations corresponds to 
a real temperature through the relation r,eai = Ttp- ImR^ks- 
Our procedure for finding a self-consistent initial condition 
thus represents an extension of the method used in Ref. [51], 
which allows us to simulate finite temperature dynamics. 

To obtain a specific initial condition, we initially set A^^ - N 
and then: (a) calculate the coefficients Ck up to the momentum 
cut-off kmax', (b) calculate A^, using Eq. (72); (c) renormalize 
the total number of atoms by reducing the condensate popu- 
lation to Nc - N - Nf. Steps (a-c) are then repeated until A^, 
converges to within 10"^ of its final value. We choose k^.^^ 
sufficiently large that Nk„.,^ < 10"^ (in addition to confirm- 
ing that the subsequent dynamics have converged as a func- 
tion of kmax)', this means that ^^ax increases substantially with 
temperature. However, for the temperatures we consider (up 
to r = 300) A;n,ax - 64 is sufficient, and the effect of the 
finite-size correction ju - is negligible. Finally, we nu- 
merically evaluate the GGPE eigenvalue Af^ for the obtained 
initial equilibrium condition. Evolution of the resulting ini- 
tial condition is accomplished using the numerical method of 
Section 111. We choose time-steps which exactly divide the di- 
mensionless kick period Tp, such that the effect of a kick can 
be accomplished via the instantaneous transformation 



UkiO) 



-IK COS 



V(0), 

Uk(0), 



vk{e)^e'''°'%k{e). 



(73) 
(74) 
(75) 
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which is applied at the exact instant of each kick. The re- 
sulting evolution is then checked for convergence in the num- 
ber of grid points, quasi-particle momentum cut-off, and time- 
step. 



C. Finite-temperature dynamics 

We now explore the dynamics of the ^-kicked rotor for a 
range of finite temperature equilibrium initial conditions. In 
addition to tracking the condensate fraction ric = Nc/N and 
non-condensate fraction n, = N,/N — variables which the 
number-conserving description gives immediate access to — 
we also introduce measures to track the presence of many- 
body, beyond-mean-field effects in the system. In Ref [51], 
the quantity 



C 



If 



(76) 



which represents a spatial average of the first-order correlation 
function g\(0,B') = (4'"'"(0')4'(6')) IN, was introduced for this 
purpose. In terms of the single-particle density matrix p{9, ff ), 
C can be defined as 



C = 



Tr [p(6»,6»')^] 



A?2 



(77) 



To an order consistent with the rest of the second-order 
number-conserving description, the single particle density 
matrix p{6, 9') is given by 

"■max 

+ 2' [(A^, + l)v,(e')vl(e)] . (78) 

The coherence measure C is thus analogous to the purity of 
a state in a single-particle system, defined by Tr[p^]/(Tr[p])^ 
(where g is the full density matrix for the single-particle sys- 
tem). Using this measure an entirely condensed state (zero 
non-condensate fraction) is akin to a pure state of a single- 
particle system, and has C = 1. The presence of non- 
condensate — either in the form of the zero-temperature quan- 
tum depletion, or in the form of thermally excited quasiparti- 
cles — guarantees that C < 1 . This is akin to a mixed state 
of a single-particle system. Physically, this can be interpreted 
as indicating the presence of high-order (many-body) corre- 
lations between the condensate and non-condensate; these are 
effectively "traced out" when one uses the number-conserving 
description to compute a single-particle density matrix, result- 
ing in p{0, ff) appearing "mixed". 

In this paper we introduce another measure which probes 
the underlying many-body correlations: the single-particle 
von Neumann entropy. This again is defined in terms of the 
single-particle density matrix, and is given by 



'sp 



-Tr 



to 
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FIG. 4. Finite-temperature evolution of the ^-kicked-rotor-BEC close 
to a nonlinear resonance explored in [50] (parameters gj = 2.5x10""*, 
A' = 10"*, Tp = 6.12, K = 0.5). Condensate and non-condensate 
fractions and n,, coherence measure C, and the single-particle von 
Neumann entropy, 5 sp are shown for initial temperatures T = 0,T = 
100, and T = 300. 



In contrast to C, which is unity for a pure condensate, 5 sp = 
for a pure condensate, and 5 sp > otherwise. Much like C, 
however, S sp gives an indication of the magnitude of the un- 
derlying many-body correlations and the degree of entangle- 
ment present in the system. It should be noted that as a direct 
measure of entanglement in a many-body system, S^p is im- 
perfect, as a proportion of its magnitude can be due to Bose- 
symmetry rather than genuine entanglement [94-96]. How- 
ever, its use as an indicator for such entanglement, as we use 
it here, is well-estabUshed [94, 97, 98]. 

To explore the finite-temperature evolution of the 5-kicked- 
rotor-BEC, we initially consider specific parameters (Tp = 
6.12, K = 0.5, gT = 2.5 X 10-\ and = 10"*) which are 
close to a nonlinear resonance; these parameters were pre- 
viously subject to detailed analysis — using the GPE, the 
first-order number-conserving description, and the second- 
order number-conserving description — in Refs. [50] and 
[51]. Figure 4 shows the evolution of the 5-kicked-rotor-BEC, 
for these parameters, resulting from three different equilib- 
rium initial conditions with temperatures T - Q, T - 100, 
and T - 300. For each initial condition the condensate frac- 
tion lie, non-condensate fraction n,, coherence measure C, and 
single-particle von Neumann entropy S^p are shown. In each 
case, the coherence measure and single-particle von Neumann 
entropy follow the non-condensate fraction very closely. This 
indicates that depletion of the condensate is the dominant fac- 
tor in creating many-body correlations and entanglement in 
this system. 
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0.80 




Dimensionless kick period, Tp 



FIG. 5. Finite-temperature shift of quantum resonances in tlie 6- 
kicked-rotor-BEC. Parameters gT = 2.5 x 10"'', A' = 2.5, and k = 0.5 
have been chosen such that T,, = 6.12 (vertical dashed line) corre- 
sponds to the nonlinear resonance explored in [50, 51] and Fig. 4. In 
particular, the finite-temperature shift in the location of this nonlinear 
resonance explains the non-resonant behavior seen at finite tempera- 
ture in Fig. 4. Interestingly, the extremely sharp cut-offs observed in 
the GPE [50] and at zero- temperature in the GGPE [5 1 ] are preserved 
at r = 300. 



The most noticeable and easily understood difference be- 
tween the temperatures is in the initial non-condensate frac- 
tion (and consequently in the coherence and entropy, since 
these generally scale with n, as explained above). This in- 
crease is to be expected, as the number of thermally excited 
quasiparticles increases sharply as a function of temperature. 
However, what is more interesting is that for increasing tem- 
peratures the impact of the periodic kicking decreases, with 
the oscillations induced by the kicks appearing to be more 
highly damped at finite temperature. This can be under- 
stood as an effect of the higher "background" thermal pop- 
ulation of the non-condensate; the large thermal population at 
higher temperatures has the effect of significantly depleting 
the condensate, and hence reducing the effective nonlinearity, 
as given by the product grN. At zero temperature, the reso- 
nant kicking period is known to shift downwards as a function 
of the product grN, which represents the effective nonlinear- 
ity [50, 51]. Hence, the presence of a thermal background, 
which effectively reduces grN, produces an upwards shift in 
the resonant kicking period. That this is the case is confirmed 
numerically in Fig. 5, in which we plot the overall response of 
the system — measured by the fractional population of all mo- 
mentum modes higher than k = Q among all atoms averaged 
over 100 kick periods (see [51]) — as a function of kicking 
period Tp and temperature. 

Another striking feature revealed by Fig. 5 is that the sharp 
cut-off at the upper limit of the shifted Talbot-time resonance 
(the large resonant area to the right of the graphs) first ob- 
served in [50] is qualitatively preserved at finite temperatures. 
This is interesting, as one might expect an effect of the in- 



creased entropy (decreased coherence) associated with finite- 
temperature to introduce more "smearing" of the resonances. 
This result is also particularly interesting in the context of ex- 
perimental realizations. For example, in an experiment with 
similar parameters to Ref [89], the unit of temperature de- 
fined above takes a value ~ 10 pK. Consequently, a future 
experiment with a temperature value in the range 1-10 nK 
would have dimensionless temperature 100 < T < 1000 in 
the units used here. Such an experiment would thus be capa- 
ble of measuring the shift in resonant frequency with temper- 
ature we predict. The preserved sharpness of the resonance 
cut-offs would help facilitate such a measurement. 

V. CONCLUSIONS 

We have described in detail a numerical method for evolv- 
ing the integro-differential equations of motion of the second- 
order number-conserving description of Gardiner and Mor- 
gan [65]. This numerical method explicitly includes the dif- 
ficult nonlinear, non-local terms which are necessary in or- 
der to preserve orthogonality between condensate and non- 
condensate, and hence gives a self-consistent treatment of 
number-dynamics and condensate - non-condensate interac- 
tions which make it ideal for studying non-equilibrium dy- 
namics in driven BEC systems at low temperatures. We have 
used this numerical method to systematically study a proto- 
typical example of such a system, the <5-kicked-rotor-BEC, at 
finite temperature. We observe that many features exhibited 
by the zero-temperature dynamics are qualitatively preserved 
at finite temperatures, and predict a shift in resonance frequen- 
cies with temperatures which could feasibly be verified by ex- 
periments using current technology. 
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Appendix A: Matrix form of the evolution operator 

1. Determination of eigenvalues 

We wish to find a general analytic expression for the matrix 
form of e"'rN(r,fM'/2/i ^ gjven position, r and time f, where 
rN(r, f) is defined by 
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( Hj^(r,t) + B(r,t) (A^i + lMi(r,/) (A^i + 2)A2(r, ?) 



rN(r,f) 



A*(r,r) 
A*(r,0 

Ai(r,0 
A2(r,0 



AM(r,t) 





















iNM + 



l)AM(r,0 iViA;(r,0 iVzA^Cr,/) 

























NMAl,ir,t)^^ 











HN(r,t) 

(Al) 



for 



and 



H^{r,t) = Vir,t) + Uo\i/f(r,t)\\ 



B(r,t) = Uoh(r,r,t)-Uo 



Nr 



A'. 



(A2) 



(A3) 



We will show that this exponentiation of FnCf, f), at a given position and time, can be accomplished in a closed analytic form for 
arbitrary M. From this point forward, we adopt the clarifying convention of omitting the expUcit position and time arguments 
r and t: all quantities appearing in our treatment should be interpreted as the complex, scalar values obtained by evaluating the 
corresponding r- and f-dependent functions at a particular position and time. 
We proceed by identifying all eigenvalues of the matrix Fn through the characteristic equation 



detCFN-^/) = 0, 

where / is the (2M + 1) x (2M + 1) identity matrix. Expanding the determinant over minors of the first row yields 



(A4) 



det(FN - ^/) = iH^+B-0 (//n -^)I-(Ni + l)Ai x 







Aj 



A* 
Al 

A2 













Hn-^ 
//n - ^ 
Hn - ^ 



+ (ATi + 1)A2 X 



A\ i/N 

A\ 








4* 
Al 

A2 



Hn-^ 
//n - ^ 
//n - ^ 























- 










(A5) 



where the choice of + and - signs in later terms is dependent on the parity of M. The determinant of the first minor is trivially 
(//gp - , and the determinant of each subsequent minor is easily computed by further expanding in minors along the row 
containing no entry of the form - ^. Carefully keeping track of signs, this yields 



.x2M-l 



det (Fn - ^/) = iH^+B-^) (//n - ^f'' + Yj^-lfiN, + 1)A,(-1)*^1a; (H^ - ^ 

M 

+ ^(-l)^"*Ar,A;(-l)^^*^iA, (//n - 



2M-1 



(A6) 
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which simplifies to 

det (Fn - m = [(Hn + B-^)(H^-0- 2] (Hn - ^f'' , (A7) 

where 



E = J](2Nk + 1)|A,|2 . (A8) 



k=-l 



Note that Eq. (A7) holds regardless of the parity of M. Consequently, 2M - 1 eigenvalues of Fn are degenerate, having value 
H^. The remaining two eigenvalues ^ solve the quadratic equation 

[HN-^f + B[//n-^]-I = 0, (A9) 

and hence are given by 

^ = //n + ^± (AlO) 



2. Determination of eigenvectors 



A set of linearly independent eigenvectors f associated with the 2M - 1 degenerate eigenvalues ^ - Hj^ must be obtained by 
solving the linear equations 



(Fn-M^^O. 

Such a set can be found by setting i/* = in Eqs. Al 1, reducing the hnear equations to 

M 

Y,[(Ni< + l)Atvl+NkAlu,]=Q. 
The required linearly independent solutions are given by 



= -6kj(Nx + l)Ai , 



where y € {2 . . . M), and 



v\ = NkAl , 

Uk = -(Jij(A^i + l)Ai , 



where j e{ \ . .. M]. 

For the remaining two eigenvalues one has to solve the linear equations 







|rN- 


V 



B IB 



(f) 



+ Z 



which can be expressed as the system 



Ak>p- 



B _ i/B 

— h 



I k=l k=l 



(All) 



(A12) 



(A13) 



(A14) 



(A15) 



(A16) 



(A17) 



(A18) 
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Rearranging the first two equations to obtain 



4* 



Uk 



(A19) 
(A20) 



provides a solution for arbitrary ifr, since substituting the above expressions into the third equation yields 

B 



and hence 



(A21) 



(A22) 



which is indeed satisfied for all lA. 

In the absence of any overriding scheme for normalizing ihe eigenvectors, a practically useful approach is to eliminate de- 
nominators in the transformation matrix P which has the eigenvectors as its columns. Choosing such a normahzation, P is given 
by: 



P = 



B+ 

Al Al (^2 + l)A2 

A* A* -(Ari + l)Ai 




(Nm + 1)Am NiAl 



A* A* 

Al Al 
A2 A2 







-{Ni + l)Ai 













<Ni + l)Ai 






iV2A; 







-(ATi + l)Ai 



[Am a 



M 



^ 

NmAI, 








-(Ni + l)Ai ) 



(A23) 



where 



(A24) 



3. Matrix Exponentiation 

In order to exponentiate the matrix Fn, one also requires the inverse, P of the transformation matrix P. Note that, since Fn 
is not Hermitian, P is not unitary, and P~^ pK However, P can nonetheless be inverted, for non-zero A^, by a lengthy series of 
elementary row operations. This yields 



P-' = 



B^-B_ 





C_(M + l)Ai C_(^2+1)A2 
C+(Ari + l)Ai C^(Ar2 + l)A2 



Ai(A^i + I)Ai 

A2(Afi+l)Ai 
W + l)Ai 



AM(Ni + \)Ai 
W + l)Ai 



A;(jV2+l)A2-£ 
(iVi + l)Ai 



A;;,(Af2 + l)A2 
W + l)Ai 

Ai(iV2+l)A2 
(Wi + DAi 

A2(jV2+l)A2 

(iVi + l)Ai 



A«(Af2 + l)A2 
(iVi + l)Ai 



C-CA^M + 1)Am C-NiA\ C-N2AI 
C^iNM + l)AM C^NiAl C^NiA*^ 



A^(AfM+l)A M 
(iVi+l)Ai 



(A'i+l)Ai 

A2(AfM+l)AM 

(iVi+l)Ai 



Am(Wm+1)Am 
(Wi + DAi 



AJiViA; A;iV2A; 



(iVi+l)Ai (Wi+ljAi 



(/Vi + l)Ai 
Ai(jVM+l)Aji, 



(/Vi + l)Ai 
AiA'iAJ-E 

(Afi+l)Ai 
AziViA 



{Ni + \)Ai 
A1N2AI 

(Ni+\)A, 
A2iV2A;-I 



(Ni+\)Ai (iVi+l)Ai 



AmWiA! 



A„A'2A! 



(iVi + l)Ai (iVi+l)Ai 



C-A^mA;, 
C+ATmA^ 

(A'i + l)Ai 



A'mNmA' 



(Ni + \)Ai 
AiNmA'^ 

(iVi + l)A, 
A2A'mA„ 

W + l)Ai 



AmNua;^- 

W + l)Ai 



(A25) 
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where 



B 



(A26) 



+ E 



Now we have Fn = PDP where D is the diagonal matrix of the eigenvalues of Fn- Consequently, we can compute the 
matrix exponential 

^-iTnStim _ p^-iDStlltip-l 



using the results we have obtained up to this point. Using the identity 

C_fi+ = -C+B = 



+ E 



and with careful treatment of all summations, one obtains 



,-!Tn5//2S . 

n _ BT^ 

COS 2 

—A*T ■ 
^1 J sin 

— A*T ■ 
^2 



—A* T ■ 
sin 



A\(Ni + l)AiT^i, 



-{N2 + i)A2r,,„ 
A*(A?2 + i)A2r„,x 



A;,(A^i + l)Air„i, A^(A^2 + 1)A2 
AiiNi + l)AiT^i, Ai(N2 + l)A2 
A2(Ni + \)AiT^, A2(N2 + 1)A2 ^mix 



-AM^sin Am(Ni + l)Air„i, Am(N2 + l)A2T„ 



-(Nm + i)AMr,i„ 

AliNM + l)AMrmix 

A^iNM + 1)AmT^. 



-n,a:t. 



\ sm 



A\NiA\T„ 
A^NiAlTn 



-N2AlT,i„ 
A\N2A*J^i, 
A*A?2A*r™ix 



AljiNM + l)AMT,^i, A*^NiAlT^,, A*^N2AIT^^ 
Ai(A?M + l)AMrmix AiN,A\T^,, AiN2A*^ 
AzCA^M + l)AMrmix A2N,A\T^,, A2N2AI 



Am{Nm + 1)AmT„ 



AMNiA\T^i^ AMN2AI 

(0 ■ 

Texp ■ 

Texp ■ 

■ 



where 



- exp 



exp 



^cos COS 



Tsin - i sin 



■ 6t 

J 

1 

2h 



exp -I 



{- 



//n + 



exp -I 



(- 



2 
B 



5t \ 
St' 

2h) 



)(l) 



+ s 



^cos 2 -^sin^ ^exp 



Thus, the action of the operator e 



-iTuSt/in 



can be succinctly expressed as 2M + 1 coupled equations: 
1 



iA(r, t + 6t/2) = \T,Ur, - ^T^^nir, t)B{r, f)j (A(r, f) - r,i„(r, f)H(r, f) , 

v*(r, f + 6tl2) = -A*(r, f)r,i„(r, f)iA(r, + re,p(r, f)v;(r, t) + T^,(r, t)E(r, t) , 
Uk{r, t + St 12) = -Ak{r, t)T,i^(r, t)if/{r, t) + T^^^ir, t)uk(r, t) + T^i^(r, f)H(r, f) , 



as used in the main text. 



(A27) 
(A28) 
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AjA^A/A^T'mix 

A2A^MA^7'niix 



A ^ M A ^ Tmix 
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A2A^MA^rmix 
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